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ABSTRACT 

Thermodynamic  and  Mechanical  Properties  of  EPON  862  with  Curing 
Agent  DETDA  by  Molecular  Simulation.  (December  2006) 

Jeremy  Lee  Tack,  B.S.,  New  Mexico  State  University 
Chair  of  Advisory  Committee:  Dr.  David  M.  Ford 

Fully  atomistic  molecular  dynamics  (MD)  simulations  were  used  to  predict  the  properties 
of  EPON  862  cross-linked  with  curing  agent  DETDA,  a  potentially  useful  epoxy  resin  for 
future  applications  of  nanocomposites.  The  properties  of  interest  were  density  (at  near¬ 
ambient  pressure  and  temperature),  glass  transition  temperature,  bulk  modulus,  and  shear 
modulus.  The  EPON  molecular  topology,  degree  of  curing,  and  MD  force-field  were 
investigated  as  variables.  The  range  of  molecular  weights  explored  was  limited  to  the 
oligomer  region,  due  to  practical  restrictions  on  model  size.  For  high  degrees  of  curing 
(greater  than  90%),  the  density  was  found  to  be  insensitive  to  the  EPON  molecular 
topology  and  precise  value  of  degree  of  curing.  Of  the  two  force-fields  that  were 
investigated,  cff91  and  COMPASS,  COMPASS  clearly  gave  more  accurate  values  for  the 
density  and  moduli  as  compared  to  experiment.  In  fact,  the  density  predicted  by 
COMPASS  was  in  excellent  agreement  with  reported  experimental  values.  However,  the 
bulk  and  shear  moduli  predicted  by  simulation  were  about  two  times  higher  than  the 
corresponding  experimental  values. 


"The  views  expressed  in  this  article  are  those  of  the  author  and  do  not  reflect  the  official 
policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense,  or  the  U.S. 
Government." 
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CHAPTER  I 
INTRODUCTION 


There  is  currently  great  interest  in  the  development  of  epoxy  resin  matrices  for  use  in  the 
composites  industry.  Private  researchers  and  the  US  military  have  particular  interest  in 
using  the  epoxy  resin  Bisphenol  F  diglycidyl  ether  (EPON  862)  with  the  curing  agent 
DETDA  for  applications  in  nanotube  composites.  When  EPON  862  is  cross-linked  with 
appropriate  curing  agents,  superior  mechanical,  adhesive,  electrical,  and  chemical 
resistance  properties  can  be  obtained  [1].  However,  creating  nanocomposites  with 
enhanced  properties  has  proven  challenging  and  there  has  been  relatively  little  modeling 
at  the  nanometer  scale  for  this  particular  system  [2], 

As  a  first  step,  we  have  carried  out  atomistic  modeling  of  neat  resins  of  EPON  862  with 
cross-linking  agent  DETDA,  with  an  emphasis  on  predicting  certain  key  thermodynamic 
and  mechanical  properties.  This  work  will  establish  a  baseline  for  future  modeling  of 
EPON-based  nanocomposites.  Molecular  topology  and  degree  of  cross-linking  are  likely 
two  critical  parameters  that  determine  the  resin’s  effectiveness  in  staying  bonded  to 
carbon  nanotubes  under  loading,  and  thus  their  effects  on  the  bulk  material  properties  are 
important  to  understand.  This  modeling  research  will  complement  ongoing  experimental 
studies  of  the  EPON-DETDA  system  using  traditional  and  nanoscale  characterization 
techniques. 


This  thesis  follows  the  style  of  Molecular  Simulation  . 
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This  main  goal  of  the  study  is  to  predict  certain  physical  and  mechanical  properties  of 
EPON  862  with  curing  agent  DETDA  using  molecular  dynamics  (MD)  simulation  and 
compare  those  predicted  with  measured  values  from  experimentation.  These  properties 
include  density,  glass  transition  temperature,  bulk  modulus,  shear  modulus,  Young’s 
modulus,  and  Poisson’s  ratio.  Another  major  issue  is  determining  which  MD  force-fields 
give  the  most  accurate  property  predictions  as  compared  to  experimental  data.  The 
sensitivity  of  the  modeling  predictions  to  assumptions  about  the  molecular  topology  and 
degree  of  curing  are  also  important  to  quantify. 
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CHAPTER  II 
BACKGROUND 


2.1  Previous  work  on  molecular  modeling  of  EPON 

Previous  literature  on  this  topic  is  very  limited.  A  combined  computational  and 
experimental  study  of  EPON  862  and  DETDA  with  single  walled  nanotubes  (SWNT) 
was  conducted  by  Zhang  et  al.  [2],  The  computational  work  employed  models  of  low 
molecular-weight  cross-linked  epoxy  species  packed  around  a  SWNT.  They  performed 
their  calculations  using  a  set  density  of  1 .2  g/cc  and  temperature  of  298  K  with  the 
COMPASS  force-field.  They  detennined  the  interfacial  bonding  energy  between  the 
epoxy  resin  and  the  nanotube  was  .1  kcal/molA2.  Using  pullout  simulations  they 
detennined  the  interfacial  shear  strength  of  up  to  75  MPa.  In  the  experimental  study  they 
found  the  uniform  dispersion  and  good  interfacial  bonding  of  the  nanotubes  in  the  epoxy 
resin  resulted  in  a  250-300%  increase  in  the  storage  modulus  with  an  addition  of  20-30 
wt%  nanotubes  [2]. 

2.2  Molecular  parameters  in  modeling  EPON-DETDA 

Before  building  the  atomistic  models,  it  is  instructive  to  consider  how  EPON  862  and  its 
DETDA-cured  polymeric  forms  are  synthesized.  In  order  to  form  EPON  862  an  epoxide 
is  combined  with  Bisphenol  F  as  seen  in  Figure  1.  Typically  the  ‘R’  group  in  the  epoxide 
is  chlorine.  Bisphenol  F  can  be  in  three  different  forms  ortho-ortho,  para-para,  and  para- 
ortho;  para-para  is  shown  in  Figure  1  [3] .  As  a  side  reaction  the  newly  formed  EPON  862 
can  also  react  with  Bisphenol  F  as  seen  in  Figure  2,  to  form  a  long  chain  with  an  epoxide 
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group  on  each  end  [4].  Since  the  side  reaction  is  driven  by  heat  it  can  be  limited.  For  this 
study  the  modified  EPON  862  in  Figure  2  is  assumed  to  have  an  ‘n’  value  of  two. 


EPON  862 


Figure  1:  Formation  of  EPON  862 


Modified  EPON  862 


OH 


Figure  2:  Addition  of  Bisphenol  F  to  EPON  862 


During  the  curing  reaction  of  EPON  862  and  DETDA  the  hydrogen  atoms  of  the  amine 
groups  of  DETDA  react  with  the  epoxide  in  EPON  862,  as  seen  in  Figure  3.  The 
resulting  molecule  can  then  react  with  up  to  three  more  epoxide  groups,  since  there  are 
still  three  more  hydrogen  atoms  left  in  the  two  amine  groups.  Figure  4  shows  the  start  of 
the  cross-linking  between  EPON  862  and  DETDA.  The  molecule  will  continue  to  react 
with  the  epoxide  and  amine  groups  expanding  in  all  directions. 
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EPON  862 


Curing  Agent  DEDTA 


EPON  862  with  Curing  Agent  DEDTA 


Figure  3:  Reaction  of  EPON  862  with  Curing  Agent  DEDTA 


EPON  862  with  Curing  Agent  DEDTA 


EPON  862 


Figure  4:  Start  of  the  Cross-linking  of  EPON  862  and  Curing  Agent 


In  summary,  when  building  the  molecular  models  of  cross-linked  EPON  862  there  are 


choices  to  be  made  about 
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(1)  Ortho/para  conformation  of  EPON  monomers 

(2)  Extent  of  side  reactions  to  produce  “extended”  EPON  monomers  (Fig  2) 

(3)  Degree  of  cross-linking  with  DETDA 

Based  on  manufacturer  data  the  ratio  of  EPON  862  to  the  curing  agent  is  100/26.4  by 
weight  [5].  When  forming  the  models  the  first  set  of  models  will  be  based  on  the 
unmodified  EPON  seen  in  Figure  1,  which  means  that  a  ratio  of  2.15  EPON  monomers  to 
one  curing  agent  will  be  used.  The  next  group  of  models  will  be  based  on  the  modified 
EPON  that  is  formed,  seen  in  Figure  2.  For  models  using  the  modified  EPON  monomer 
the  ratio  of  EPON  to  curing  agent  will  decrease  since  the  molecular  weight  of  the 
modified  EPON  is  higher,  therefore  a  ratio  of  2  to  1  will  be  used. 

2.3  Molecular  modeling  techniques 
2.3.1  Software 

There  are  several  software  packages  that  are  available  for  molecular  simulation  that  could 
be  used.  The  two  that  will  be  used  in  the  current  work  are  LAMMPS  and  Materials 
Studio.  LAMMPS  is  a  classical  molecular  dynamics  (MD)  code  that  can  model  systems 
in  either  a  liquid,  gaseous,  or  solid  state.  The  code  was  written  by  Steve  Plimpton  and 
was  developed  at  Sandia  National  Laboratories  with  funding  from  the  Department  of 
Energy.  It  is  an  open-source  code  that  is  freely  available  [6]. 

Materials  Studio  was  developed  by  Accelrys  Software  Inc.  The  three  main  products  of 
Materials  Studio  that  will  be  used  are  Amorphous  Builder,  Discover,  and  COMPASS. 
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Amorphous  Builder  is  used  for  building  periodic  cells  of  systems  and  designed  especially 
for  liquids  and  amorphous  polymers.  Discover  is  the  MD  simulation  tool  used  in 
Materials  Studio.  COMPASS  stands  for  Condensed-phase  Optimized  Molecular 
Potentials  for  Atomistic  Simulation  Studies.  It  is  a  force-field  that  is  only  available  on 
Accelrys  software  [7]. 

2.3.2  Molecular  dynamics 

In  classical  molecular  simulation,  a  model  system  is  built  at  the  atomistic  level  with 
prescribed  potentials  (the  force-field)  acting  between  the  atoms.  These  interactions 
consist  of  site-site  interactions,  such  as  van  der  Waals  dispersion  and  Coulombic  forces, 
and  intramolecular  forces,  such  as  chemical  bonds,  angle  bending  and  dihedral  torsional 
barriers.  Models  typically  span  on  the  order  of  1  to  10  mn  on  a  side,  with  periodic 
boundary  conditions  in  one  or  more  of  the  dimensions  [8].  The  term  periodic  boundary 
condition  refers  to  the  simulation  of  structures  consisting  of  a  periodic  lattice  of  identical 
subunits.  During  the  course  of  the  simulation  as  a  molecule  moves  in  the  original  cell,  its 
image  in  the  other  cells  moves  in  exactly  the  same  way [8].  Molecular  and  macroscopic 
properties  of  interest  can  be  averaged  over  a  molecular  dynamics  (MD)  trajectory,  which 
is  simply  an  integration  of  Newton’s  equations  of  motion  for  all  of  the  atoms  in  the 
system.  Time  spans  on  the  order  of  nanoseconds  may  be  routinely  accessed  with  desktop 
computers  [8,9]. 

An  important  issue  in  molecular  dynamics  simulations  is  the  summation  of  the  pairwise 
(van  der  Waals  and  Coulombic)  potential  functions.  In  principle,  computing  their  force 
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2 

contribution  to  the  equations  of  motion  requires  (N  -N)/2  tenns,  where  N  is  the  total 
number  of  atoms.  For  this  study  around  10,000  atoms  are  used  for  each  simulation, 
which  equates  to  approximately  50  million  terms.  Typically  the  problem  is  made 
tractable  by  splitting  these  forces  into  short-range  and  long-range  contributions.  A  cut¬ 
off  radius,  rcut,  is  introduced  which  represents  an  effective  maximum  range  for  the  van 
der  Waals  forces.  For  distances  that  are  smaller  than  rcut,  the  van  der  Waals  forces  are 
calculated  using  simple  summation;  the  limited  range  significantly  reduces  the  number  of 
pairs  that  must  be  considered.  For  the  long-range  (Coulombic)  contributions,  methods 
such  as  Ewald  summation  or  particle -particle  and  particle-mesh  (PPPM)  are  the  common 
methods  [10].  Ewald  summation  decomposes  the  Coulombic  forces  into  real  space  and 
Fourier  space  contributions  using  the  periodicity  imposed  by  the  boundary  conditions. 

In  order  to  extract  infonnation  from  molecular  simulations  a  convenient  ensemble  should 
be  chosen.  The  two  that  are  used  in  this  study  are  NPT  and  NVT.  NPT  holds  the  number 
of  atoms,  pressure,  and  temperature  of  the  system  constant.  The  NVT  hold  the  number  of 
atoms,  volume,  and  temperature  of  the  system  constant.  The  NPT  ensemble  is  the  most 
useful  in  this  study  since  the  mechanical  properties  are  known,  from  experiments,  at 
certain  temperature  and  pressure  values. 

2.4  Temperature  and  pressure  control  methods 
2.4.1  Temperature  control 

For  controlling  the  temperature  in  a  MD  simulation  there  are  four  main  methods  velocity 
scale,  Nose,  Andersen,  and  Berendsen.  Velocity  scale  controls  the  kinetic  temperature  of 
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a  system  and  brings  it  to  equilibrium  by  maintaining  the  temperature  within  a  given  range 
of  the  target  temperature.  When  the  temperature  drifts  outside  of  the  specified  range  the 
atomic  velocities  are  simply  scaled  back  into  range.  Velocity  rescaling  method  is 
physically  unrealistic  and  should  not  be  used  for  true  temperature  control.  For  Nose, 
Andersen,  and  Berendsen  control  the  thermodynamic  temperature  and  generate  the 
correct  statistical  ensemble,  by  allowing  the  system  to  exchange  energy  with  a  heat  bath 

m. 


2.4.2  Pressure  control 

For  controlling  the  pressure  in  a  MD  simulation  there  are  four  main  methods  Parinello, 
Andersen,  Nose  and  Berendsen.  The  Parinello  method  is  useful  for  studying  the  stress- 
strain  relationship  in  materials.  Both  the  shape  and  the  volume  of  the  cell  can  change, 
which  allows  the  internal  stress  of  the  system  to  match  the  externally  applied  stress.  The 
Andersen  and  Nose  method  allows  changing  of  the  volume  of  the  cell  but  keeps  the  shape 
fixed  by  allowing  the  cell  to  change  isotropically.  The  Berendsen  method  involves 
changing  the  pressure  by  altering  the  coordinates  of  the  particles  and  the  size  of  the  cell, 
the  shape  is  also  kept  fixed  [7]. 

2.5  Properties  predicted 
2.5.1.  Density 

Density  at  fixed  temperature  and  pressure  is  a  key  output  being  measured.  The  density  is 
obtained  from  an  NPT  simulation  by  dividing  the  total  mass  of  atoms  in  the  simulation 
cell  by  the  average  cell  volume  over  the  MD  run. 
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2.5.2  Bulk  and  shear  moduli 

The  design  of  this  portion  of  the  study  begins  with  the  classical  equations  for  the  shear 
modulus  (G)  and  bulk  modulus  ( B )  in  terms  of  the  first  and  second  Lame  constants  X  and 
p,  as  shown  in  Equations  1  and  2  respectively  [11].  The  shear  modulus,  also  known  as 
the  modulus  of  rigidity,  is  the  measure  of  elastic  deformation  of  a  body  in  which  an 
applied  force  determines  the  shape  of  the  body  [12].  The  bulk  modulus  is  the  ratio  of  the 
change  in  pressure  acting  on  a  volume  to  the  fractional  change  in  volume  [12]. 

G  =  ju  Eq  1 

2 

B  =  A  +  —/i  Eq2 

The  simulations  run  in  this  study  seek  to  find  B  and  G  in  Materials  Studio  and  use  those 
values  to  find  the  Lame  constants  X  and  p.  The  bulk  modulus  can  be  easily  determined 
from  Equation  3  by  simply  introducing  a  high,  e.g.  5000  atm,  hydrostatic  pressure  to  the 
NPT  simulation  and  observing  the  new  average  cell  volume  upon  equilibration  at  this 
new  pressure. 


B  = 


Eq  3 


Applying  hydrostatic  pressure  induces  no  shear  stress  in  the  control  volume,  so  it 
provides  no  infonnation  on  G.  Therefore,  at  this  point,  additional  NPT  simulation  runs 
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with  an  induced  XZ  shear  stress  were  done  with  each  of  the  samples  to  determine  G.  A 
free  body  diagram  of  the  loading  configuration  appears  in  Figure  5. 


N 


X 


Figure  5:  FBD  of  Loading  Configuration 


The  method  for  detennining  the  shear  modulus  from  the  data  in  Materials  Studio  takes  a 
known  value  of  shear  stress  as  well  as  initial  lattice  parameters  (a,b,c,a,P,y)  and  looks  at 
how  the  lattice  deforms  over  time  to  find  the  shear  modulus.  A  method  for  doing  this 
calculation  in  free  space  (xyz)  coordinates  exists  but,  since  the  lattice  parameter 
directions  do  not  necessarily  correspond  with  the  x,  y,  and  z  axes,  they  must  be  arbitrarily 
aligned.  This  alignment  is  shown  below  in  Figure  6  where  the  c  direction  sits  aligned 
with  the  z  direction,  the  b  direction  is  rotated  from  the  z  axis  by  the  lattice  parameter  a, 
and  the  direction  is  aligned  off  of  the  z  axis  by  the  angle  P  and  the  other  two  directions 


with  a  new  angle  y  . 
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Figure  6:  Alignment  of  Lattice  Parameters 

With  this  alignment  in  place,  it  becomes  possible  to  convert  the  lattice  parameters  a,  b, 
and  c  into  the  Cartesian  coordinates  x,  y,  and  z  with  the  equations  in  Equation  4a  and  4b. 


ax  =  a*  sin  /?  *  sin  y  ay  =  a  *  sin  /?  *  cos  y  az  =  a*  cos  /?  bx  =  0  c  =  0 
by=b*sin(a)  bz  =  b*cos(a)  c=c  c=Q 


Eq  4a 


-  cos  a  *  cos/?  +  cosy 

cos  y  = - - - 

sin  a*  sin  /?  gq  4^ 

sin y  =  -y/l  -cos2  y 


These  Cartesian  parameters  determine  the  lattice  matrix  H  (Equation  5)  which  can  be 
used  to  directly  determine  the  3D  strain  matrix  from  Equation  6  [13].  Shear  strain  in  the 
xz  direction  can  be  directly  taken  from  the  strain  matrix  and  substituted  into  Equation  7 
along  with  the  known  shear  stress  to  find  the  shear  modulus. 
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ax 

bx 

Cx 

ay 

by 

Cy 

a: 

bx 

C: 

Eq  5 


Eq  6 


G  = 


£ 


xz 


Eq  7 


Where  Ho  denotes  the  converted  lattice  parameter  matrix  at  the  initial  simulation 
conditions  [14]. 

2.5.3  Glass  transition  temperature 

The  method  for  detennining  the  glass  transition  temperature  (Tg)  uses  density  or  volume 
data  from  different  temperature  values  over  a  given  range.  Temperature  is  then  plotted 
versus  density.  Before  and  after  the  Tg  the  slope  of  the  plot  should  be  a  constant  value. 
This  allows  the  Tg  to  be  detennined  by  finding  the  point  where  the  slope  changes. 

2.6  Run  details 

2.6.1  Building  the  atomistic  models 

Once  the  configuration  of  the  EPON  monomer;  ortho-ortho,  para-para,  or  para-ortho;  and 
whether  to  use  either  the  nonnal  EPON  monomer,  Figure  1,  or  the  modified  EPON 
monomer,  Figure  2,  is  chosen  as  described  in  section  2.2,  a  model  is  built  comprising  the 
desired  number  of  each  type  of  molecule  packed  into  a  simulation  volume  with  periodic 
boundary  conditions.  For  this  purpose  Materials  Studio’s  Amorphous  Builder  package 
was  used.  Amorphous  Builder  inserts  the  molecules  that  you  select  into  a  periodic  cell 
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with  a  density  that  you  specify,  usually  lower  than  the  target  density  for  the  system. 

Once  the  periodic  cell  is  made  the  density  can  be  adjusted  by  changing  the  periodic  cell 
parameters  in  small  increments,  minimizing  each  time  which  prevents  overlaps  between 
atoms,  until  the  target  density  is  reached  making  sure  that  cell  is  kept  in  a  box 
configuration.  The  system  is  then  minimized,  by  optimizing  the  molecular  structure, 
until  the  system  total  energy  is  at  a  minimum.  Now  the  system  is  ready  for  dynamics  to 
be  run. 

2.6.2.  Equilibration  and  production  runs 

Now  the  temperature  and  pressure  need  to  be  stabilized  before  the  production  runs  can 
start.  Using  a  NVT  ensemble  and  Berendsen  temperature  control  the  system  is  run  for  at 
least  20ps  to  stabilize  the  temperature  at  298K.  Next  a  NPT  ensemble  using  Berendsen 
for  both  temperature  and  pressure  control  the  system  is  run  for  another  20ps  to  stabilize 
the  pressure  at  latm.  Now  the  system  is  ready  for  production  runs. 

For  production  runs  either  Materials  Studio  or  LAMMPS  was  used.  See  Table  1  for 
specific  information  on  force-field  and  control  methods  used.  For  both  programs  Ewald 
summation  will  be  used  for  long-range  Coulombic  forces  and  Van  der  Waal  summation 
for  short-range.  See  Table  2  for  specific  information  on  the  adjustable  parameters  used. 
The  simulations  will  be  run  for  500  ps  to  start  with,  until  it  is  detennined  when 
equilibrium  is  reached,  after  which  the  total  time  run  will  be  reduced,  still  above  the  time 
to  reach  equilibrium,  in  order  to  save  computational  time. 
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Starting  the  runs  is  different  for  each  program  that  is  used.  In  Materials  Studio,  since  the 
models  where  built  in  this  program  the  only  thing  needed  is  to  start  Discover  to  run  them. 
For  LAMMPS  Discover  is  still  used,  but  instead  of  running  the  simulation  it  needs  to  be 
saved  to  a  file,  making  sure  the  correct  force-field  is  selected.  Then  making  use  of  a 
converter  that  comes  with  LAMMPS  the  file  is  put  into  the  correct  fonnat  for  use.  Since 
LAMMPS  is  a  UNIX  program,  a  script  file  is  used  to  start  the  simulation. 


Table  1:  Control  Methods  for  Software  Used 


Materials  Studio 

LAMMPS 

Force-field 

COMPASS 

Cff91 

Temperature  Control 

Berendsen 

Nose/Hoover 

Pressure  Control 

Berendsen 

Nose  Hoover 

Short  Range  Summation 

van  der  Waals 

van  der  Waals 

Long  Range  Coulombic 

Ewald 

Ewald 

Table  2:  Adjustab 

e  Parameters  Used  for  Physical  and  Mechanical  Properties 

Temperature  (K) 

Pressure  (atm) 

Shear  Stress  (atm) 

Density 

298 

1 

0 

Bulk  Modulus 

298 

5000 

0 

Shear  Modulus 

298 

1 

100 
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CHAPTER  III 
RESULTS 


3.1  Macromolecular  topology/degree  of  curing  effects  on  density 

Starting  with  the  assumption  that  just  normal  EPON  862  monomers  and  curing  agent 
DETDA  are  used,  without  the  modified  EPON  862  structure  seen  in  Figure  2,  several 
different  atomistic  models  were  built  and  used  with  MD  to  predict  the  density  at  298  K 
and  1  atm. 


Figure  7:  EPON  862  and  Curing  Agent  Cross-Linked  in  a  Square  Configuration 


Figure  8:  EPON  862  and  Curing  Agent  Cross-Linked  in  a  Triangle  Configuration 
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Figure  9:  Two  Views  of  Square  Model  System  Studied,  (a)  Represents  the  default  view 
where  chemically  bonded  species  are  represented  in  extended  form,  (b)  The  same  model 
that  represents  the  periodic  boundary  conditions  implemented  where  the  atoms  from 
neighboring  cells  are  displayed. 


The  square,  Figure  7,  and  triangle,  Figure  8,  runs  both  have  10  molecules,  each  composed 
of  13  EPON  monomers  and  6  curing  agents  configured  in  a  square  or  triangle  shape 
respectively.  Figure  9  shows  two  views  of  the  square  model  system,  built  with  the 
Amorphous  Builder  software  as  described  in  section  2.6.1,  where  the  one  cell 
representation,  (b),  shows  an  almost  uniform  distribution  of  atoms  filling  the  simulation 
cell  of  an  amorphous  system.  Since  the  individual  molecules  (Figures  7  and  8)  are 
somewhat  two-dimensional  by  construction,  one  concern  is  that  the  final  packed  models 
(Figure  9)  might  have  unrealistic  structural  anisotropies  (although  the  techniques  of 
Amorphous  Builder  are  designed  to  maximize  the  isotropy).  From  visual  inspection  of 
the  packed  models,  using  representations  like  Figure  9  and  also  the  free-volume 
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visualization  tool  in  Materials  Studio,  there  was  no  obvious  indication  of  anisotropic 
conditions.  However,  this  does  not  guarantee  that  some  level  of  anisotropy  is  not  present. 

The  partially  cured  run  has  9  molecules,  each  composed  of  13  EPON  monomers  and  6 
curing  agents  in  a  square  configuration  and  13  EPON  monomers  and  6  curing  agent 
molecules  that  are  not  bonded  to  each  other,  which  corresponded  to  a  90  percent  cured 
sample.  The  liquid  run  was  a  set  of  EPON  monomers  and  curing  agents  that  are  not 
bonded  to  each  other.  Table  3  summarizes  the  data. 


Table  3:  Comparison  of  Density  1 

for  Different  Model  Configurations 

cured90 

liquid 

square 

triangle 

Density  Range 

1.061-1.078 

1.002-1.024 

1.036-1.061 

1.018-1.049 

Average 

1.067 

1.011 

1.054 

1.039 

Standard  Deviation 

0.0088 

0.0085 

0.012 

0.014 

Time  (ps) 

600 

1000 

600 

900 

Number  of  runs 

5 

5 

4 

4 

These  four  models  show  that  the  molecular  topology  and  degree  of  curing  (above  90%) 
does  not  have  a  statistically  significant  effect  on  the  density.  On  the  other  hand,  a 
completely  unbonded,  (liquid)  model  does  predict  a  significant  lower  density.  Since  the 
molecular  topology  does  not  have  a  significant  effect,  the  square  configuration  will  be 
used  for  the  rest  of  the  tests.  All  of  these  run  also  showed  that  they  where  at  equilibrium 
before  75  ps,  therefore  to  reduce  computer  time  the  rest  of  the  runs  will  be  run  for  100  ps. 
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Since  it  is  known  that  EPON  862  can  have  a  modified  structure  as  seen  in  Figure  2, 
another  set  of  atomistic  models  were  built  using  the  modified  structure.  All  runs  were 
obtained  with  LAMMPS  using  the  force-field  cff91  at  298K  and  1  atm. 

The  modified  runs  have  10  molecules,  each  composed  of  7  normal  EPON  monomers  and 
1  modified  EPON  monomers  and  4  curing  agents  configured  in  a  square.  The  modified 
cured  run  has  9  molecules,  each  composed  of  7  normal  EPON  monomers  and  1  modified 
EPON  monomers  and  4  curing  agents  configured  in  a  square  and  7  normal  EPON 
monomers  and  1  modified  EPON  monomer  and  4  curing  agent  molecules  that  are  not 
bonded  to  each  other.  Half  of  the  models  used  the  ortho-ortho  configuration  for  EPON 
862  and  the  other  used  the  para-para  configuration.  Experimental  data  from  manufacture 
web  site  [5],  Table  4  summarizes  the  data. 


Changing  to  the  modified  EPON  monomer  led  to  only  a  two  percent  increase  in  the 
density.  Whereas  changing  the  configuration  from  ortho-ortho  to  para-para  did  not 
significantly  change  the  range  of  density  values. 


Table  4:  Comparison  of  Different  Configurations  of  EPON  862 


Ortho-Ortho 

Para-Para 

Cured 

Experimental 

Density  Range 

1.054-1.082 

1.058-1.083 

1.068-1.074 

1.17-1.2 

Average 

1.073 

1.075 

1.07 

Standard  Deviation 

0.0114 

0.0105 

0.0028 

Time  (ps) 

100 

100 

100 

Number  of  runs 

5 

5 

3 
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3.2  Force-field  effects  on  density 

The  quality  of  the  force-field  that  is  applied  in  MD  simulations  is  critical  to  the  accuracy 
of  the  predications.  In  addition  to  the  cff91  force-field,  the  COMPASS  force-field  was 
also  evaluated.  The  force-field  cff91  is  freely  available  and  transferable  to  any  MD 
simulation  software,  unlike  COMPASS  which  is  proprietary  to  Accelrys  Inc  and  can  only 
be  run  on  their  software  (e.g.  Materials  Studio  or  Cerius2). 

Due  to  the  fact  that  Materials  Studio  is  run  on  a  PC  environment,  this  limits  the  size  of 
the  model  and  how  fast  it  runs.  All  runs  on  Materials  Studio  were  reduced  from  10 
molecules  to  5,  but  otherwise  where  the  same  as  the  modified  runs  using  LAMMPS. 
Experimental  data  are  from  the  manufacturer’s  web  site  [5].  Table  5  summarizes  the 
data. 


Table  5:  Comparison  of  Density  for  Different  Force-fields 


Cff91 

COMPASS 

Experimental 

Density  Range 

1.054-1.082 

1.123-1.1345 

1.17-1.2 

Average 

1.074 

1.13 

Standard  Deviation 

0.0104 

0.0061 

From  Table  4  it  is  clear  that  by  changing  the  force-field  the  density  has  increased  by  over 
five  percent.  This  shows  that  the  force-field  has  a  significant  effect  on  the  model  that  is 
being  used.  The  density  is  now  3-6%  lower  than  the  experimental  values  using 


COMPASS  versus  8-11%  lower  with  cff91. 
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3.3  Mechanical  properties 

First,  the  bulk  modulus  test  was  run  using  an  NPT  ensemble  with  temperature  set  to  298K 
and  pressure  set  to  equilibrate  at  5000  atm.  The  equilibrating  pressure  changes  the 
volume  of  the  periodic  cell  by  an  amount  equal  to  AV  in  Equation  3.  Knowing  how  the 
pressure  changed  and  the  initial  and  final  dimensions  of  the  control  volume,  the  bulk 
modulus  can  be  calculated.  Experimental  bulk  modulus  data  is  from  the  research  group 
of  D.  Lagoudas  [15].  Table  6  summarizes  the  results. 


Table  6:  Comparison  of  Bulk  Modulus  for  Different  Force-fields 


Cff91 

COMPASS 

Experimental 

Bulk  Modulus(GPa)  Range 

5.64-6.036 

4.84-5.54 

Average 

5.79 

5.18 

2.9 

Standard  Deviation 

0.213 

0.535 

It  is  clear  that  the  COMPASS  force-field  estimates  that  bulk  modulus  better  than  cff91 
does,  but  both  are  clearly  higher,  1.5-2x,  than  the  experimental  value.  This  is  because  the 
simulation  system  represents  an  ideal  system,  whereas  the  experimental  system  may  have 
microscopic  or  even  macroscopic  defects.  So  the  simulated  values  can  be  seen  as  the 
upper  limit  on  what  the  value  can  be  [16]. 

Next,  the  shear  test  was  run  using  another  NPT  ensemble  with  temperature  of  298K. 
Pressure  control  is  the  key  factor  is  this  run.  The  pressure  was  set  to  one  atm,  but  a  shear 
stress  of  100  atm  was  added  in  the  XZ  cell  direction  and  Parinello  barostat  control  was 


used.  The  change  to  a  Parinello  control  system  for  pressure  was  due  to  the  fact  that 
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Berendsen  control  keeps  the  lattice  in  the  original  cubic  shape,  but  to  calculate  the  shear 
modulus  the  lattice  shape  needs  to  be  able  to  change  with  the  shear.  The  ensemble 
provides  the  lattice  parameter  values  necessary  in  the  calculation  described  in  section 
2.5.2.  Occasionally  when  running  the  system  with  a  shear  stress  it  was  noted  that 
erroneous  results  were  received.  On  average  the  erroneous  results  occurred  in  about  ten 
percent  of  the  runs  conducted.  The  usual  result  was  that  the  system  would  not  run  to 
completion,  i.e.  the  system  would  crash.  Occasionally  the  results  would  be  an  unrealistic 
result,  usually  extremely  high  or  a  negative  number.  The  extremely  high  results  usually 
were  caused  when  the  system  would  be  on  the  verge  of  crashing,  but  the  run  finished 
before  this  happened  causing  the  values  reported  to  be  very  high. 

With  the  values  of  B  and  G  known,  calculation  of  the  Lame  constants  as  well  as 
Poisson’s  ratio(v)  and  Young’s  modulus  (E)  became  possible  using  Equations  1  and  2 
above  as  well  as  the  two  equations  below  in  Equation  8[1 1],  Poisson's  ratio  is  the  ratio  of 
transverse  contraction  strain  to  longitudinal  extension  strain  in  the  direction  of  stretching 
force  [17].  Young’s  modulus  is  a  measure  of  the  stiffness  of  a  given  material.  It  is 
defined  as  the  ratio,  for  small  strains,  of  the  rate  of  change  of  stress  with  strain  [18], 


E  = 


M 


3  A  +  2  fi 
A  +  /u 


A 

1  ~  2(A  +  ju) 


Eq  8 


From  experiments  in  the  research  group  of  D.  Lagoudas,  the  bulk,  storage  (E’),  and  loss 
(E”)  moduli  are  known  [15]  see  Table  7. 
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Table  7:  Experimental  Values  of  Storage  and  Loss  ] 

VIodulus 

Experimental 

Calculated 

Storage  (E’) 

Loss  (E") 

Young’s  (E) 

Sample  1 

2.843 

0.21 

2.85 

Sample  2 

2.913 

0.223 

2.92 

Sample  3 

1.908 

0.113 

1.91 

Average 

2.555 

0.182 

2.59 

Standard  Deviation 

0.561 

0.0601 

0.564 

The  Young’s  modulus  can  be  calculated  from  the  storage  and  loss  modulus,  see  Table  7, 
by  using  Equation  9  [19].  Having  the  Bulk  and  Young’s  modulus  allows  the  calculation 
of  the  shear  modulus  and  Poisson’s  ratio. 

E  =  E'+iE' '  |£|  =  J(E')2  +  {E")1  Eq  9 

Final  calculated  values  are  presented  below  in  Table  8. 


Table  8:  Comparison  oi 

f  Mechanical  Properties 

k(GPa) 

p(GPa) 

Bulk(GPa) 

Shear(GPa) 

Young ’s(GPa) 

Poisson 

COMPASS 

4.81 

.557 

5.18 

.557 

1.61 

.448 

Experimental 

2.27 

.946 

2.9 

.946 

2.56 

.353 
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The  reason  that  only  COMPASS  results  are  shown  is  that  LAMMPS  is  unable  to  predict 
the  shear  modulus  using  the  method  described  above. 

3.4  Glass  transition  temperature 

The  glass  transition  temperature  (Tg)  can  be  obtained  by  generating  density  data  at 
different  temperatures  and  looking  for  a  change  in  the  slope  of  a  density  vs.  temperature 
plot  [20,21].  See  Figure  10 


Glass  Transition  Temperature 


♦  COMPASS 
■  Cff91 


250  350  450  550 

Temperature  (K) 

Figure  10:  Plot  of  Density  vs.  Temperature 

By  looking  at  Figure  10  the  Tg  can  be  seen  to  be  between  350  and  400K  for  COMPASS 
and  between  400  and  425K  for  cff9 1 .  Experimentally  the  value  of  Tg  have  been 
reported,  by  the  research  groups  of  D.  Lagoudas  and  FI.  J.  Sue,  between  418  and  434K 
[15,22].  From  this  it  seems  that  for  estimating  the  Tg,  the  cff9  I  force-field  is  the  better 


choice. 
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CHAPTER  IV 
CONCLUSIONS 


This  study  of  EPON  862  with  curing  agent  DETDA  using  molecular  dynamics  simulation 
focused  on  three  main  points:  predict  key  physical  and  mechanical  properties;  evaluate 
the  effects  of  molecular  topology  on  the  predictions;  and  detennine  which  type  of  force- 
field  gives  the  most  realistic  results  as  compared  to  experiment. 

The  physical  and  mechanical  properties  that  were  predicted  and  compared  with 
experimental  values  were  density,  glass  transition  temperature,  bulk  modulus,  shear 
modulus,  Young’s  modulus,  and  Poisson’s  ratio.  In  the  best  cases,  the  predicted  density 
was  within  6%  of  the  experimental  value  and  the  predicted  Tg  was  within  3%.  The 
modulus  is  where  molecular  dynamics  has  a  problem  predicting  the  values.  The 
calculated  bulk  modulus  was  at  least  1.5  times  higher  than  the  experimental  value.  This 
is  because  the  simulated  system  represents  an  ideal  system,  whereas  the  experimental 
may  have  microscopic  or  even  macroscopic  defects.  The  Young’s  modulus  was  about 
0.6  times  the  experimental  value  and  Poisson’s  ratio  was  about  1.3  times  higher  than  the 
experimental  value.  The  reasons  that  the  Young’s  modulus  is  lower  than  the 
experimental  value  are  likely  (1)  slightly  non-isotropic  systems  and  (2)  the  models  being 
used  are  smaller  than  those  used  experimentally. 

For  the  molecular  topology  effects,  the  ortho/para  confonnation,  degree  of  cross-linking 
with  DETDA,  and  the  extent  of  side  reactions  were  examined.  EPON  862  can  be  in 
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ortho-ortho,  para-para,  or  para-ortho  types  of  conformation.  It  was  found  that  the  type  of 
conformation  had  no  statistical  impact  on  the  predicted  density  of  EPON  at  near-ambient 
conditions.  The  amount  of  cross-linking  directly  impacts  the  physical  properties. 

Samples  with  no  cross-linking,  where  none  of  the  molecules  are  bonded  to  each  other, 
were  compared  to  oligomeric,  i.e.  a  polymer  with  only  a  few  monomers,  samples  with 
cross-linking  based  on  manufacturer  data;  a  5-6%  increase  in  the  density  with  cross- 
linking  was  observed  using  the  cff9 1  force-field.  For  the  cross-linked  samples,  the 
particular  configuration  of  the  oligomers  had  little  impact  on  the  predicted  density.  The 
side  reaction,  which  added  an  extra  Bisphenol  F  group  to  EPON  862,  showed  only  a 
slight  increase  in  the  density,  about  2%.  Overall  the  aspect  of  molecular  topology  that 
most  affected  the  physical  properties  was  the  amount  of  cross-linking. 

Two  force-fields  were  evaluated  during  this  study:  cff91  using  EAMMPS,  and 
COMPASS  using  Materials  Studio.  The  cff9 1  force-field  was  found  to  be  more  effective 
only  in  calculating  the  glass  transition  temperature.  For  everything  else  COMPASS  was 
more  accurate,  as  judged  by  comparison  with  experimental  data.  The  drawback  of  using 
COMPASS  is  that  it  is  proprietary  to  Accelrys,  whereas  cff91  can  be  used  with  any 
software  package. 
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CHAPTER  V 
FUTURE  WORK 


Since  the  mechanical  properties  that  are  being  calculated  are  not  matching  what  is  being 
found  experimental  something  needs  to  be  done.  Due  to  the  fact  that  the  force-field 
employed  is  critical  to  the  accuracy  of  the  predictions  of  the  simulations,  better  force- 
fields  need  to  be  found.  One  thing  that  can  be  done  is  to  develop  a  new  force-field  from 
scratch  that  predicts  the  mechanical  properties  more  accurately.  Additionally  existing 
force-fields  could  be  modified.  While  this  should  bring  the  mechanical  properties  closer 
to  the  experimental  values  it  will  not  solve  all  the  problems. 

Microscopic  and/or  macroscopic  defect  could  also  be  affecting  the  mechanical  properties. 
Microscopic  defects  could  be  included  in  the  current  models  that  are  being  used,  but  not 
knowing  what  type  of  defects,  and  how  many,  that  could  be  found  in  this  system  limit  the 
testing  capabilities.  Macroscopic  defects  cannot  be  tested  in  our  models  since  the  defect 
would  be  as  large  as  the  total  model.  Macroscopic  defects  would  need  to  be  tested  using 
mesoscale  dynamics.  A  way  to  bridge  molecular  and  mesoscale  dynamics  needs  to  found 
in  order  for  the  full  effects  of  defects  to  be  tested. 

Now  that  the  structural  behavior  and  the  mechanical  properties  are  known,  the  next  step 
in  this  project  is  to  start  adding  nano  tubes  into  the  molecular  simulation.  The  nano  tube 
should  be  surrounded  by  cross-linked  EPON  862  with  DETDA.  Several  problems  arise 
when  doing  this.  First  is  since  Amorphous  Builder  in  Materials  Studio  inserts  molecules 
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randomly  into  a  periodic  cell  the  nanotube  might  not  be  where  it  should  be.  Therefore 
several  iterations  will  be  necessary  until  the  nanotube  is  positioned  correctly.  Once  the 
nanotube  is  in  the  correct  position,  nonnal  molecular  dynamics  will  not  be  a  problem 
until  high  pressure  is  added  to  the  system.  Due  to  the  cylindrical  structure  of  nanotubes 
problems  arise  when  adding  high  pressure  to  the  system,  for  calculation  of  the  Bulk 
modulus.  When  the  nano  tube  is  compressed  it  tends  loss  its  cylindrical  shape,  when 
using  the  COMPASS  force-field.  Other  force-fields  might  not  have  this  same  problem 
and  needs  to  be  looked  into.  A  method  will  need  to  be  developed,  using  current  force- 
fields  to  make  sure  that  the  nanotubes  retain  their  shape  but  still  allow  for  heavy 
compression  of  around  5000  atm.  Adding  a  shear  stress  to  the  system  might  cause  the 
same  type  of  problem  and  will  need  to  be  looked  into. 

Once  the  solutions  to  these  problems  have  been  solved  additional  nanotubes  will  be 
added  one  at  a  time  to  the  system  until  they  no  longer  give  a  significant  increase  to  the 
mechanical  and  structural  properties  of  the  system.  Knowing  the  upper  limit  on  the 
number  of  nanotubes  that  can  be  added  to  the  system  will  allow  cost  versus  impact 
projections  to  be  made.  Once  the  optimum  number  of  nano  tubes  to  be  added  to  the 
system  is  known,  it  will  give  the  experimentalists  a  good  starting  point  for  their  testing. 
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